Hydrodynamic excitations of Bose Condensates in anisotropic traps 
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The collective excitations of Bose condensates in anisotropic axially symmetric harmonic traps are 
investigated in the hydrodynamic and Thomas-Fermi limit. We identify an additional conserved 
quantity, besides the axial angular momentum and the total energy, and separate the wave equation 
in elliptic coordinates. The solution is thereby reduced to the algebraic problem of diagonalizing 
finite dimensional matrices. The classical quasi-particle dynamics in the local density approximation 
for energies of the order of the chemical potential is shown to be chaotic. 
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The new Bose condensates of alkali atoms in magnetic 
traps offer a unique way to investigate the low-lying 
collective excitations in Bose condensates [^|dT6|. Exper- 
imentally collective modes with a given symmetry have 
been excited by time-dependent modulations of the trap- 
ping potential, and their evolution has been followed in 
real time by measurements of the resulting shape oscilla- 
tions of the condensates. The measurements performed 
so far have involved turning off the trap after a given time 
j§-||, but in the future they could even be performed 
non-destructively by elastic off-resonant light scattering 
0. Theoretically the collective modes have been ana- 
lyzed by using the Bogoliubov-equations or by lineariz- 
ing the time-dependent Gross-Pitaevskii equation around 
the time- independent condensate and solving these equa- 
tions numerically |^,^| or analytically in various approx- 
imations pTcj— [T^] . Very good agreement between the nu- 
merical and the experimental results has been found. 

In a seminal paper Stringari |l6| has shown how the 
coupled wave equations for the collective excitations are 
simplified in the hydrodynamic limit to become a sin- 
gle second-order wave equation for density waves, and he 
obtained analytical solutions for all its modes in spheri- 
cally symmetric harmonic traps, and, remarkably, also 
for some of its modes in axially symmetric harmonic 
traps. The latter are particularly important, because all 
experiments have been performed with traps of this sym- 
metry |~§. 

In the present paper it is our goal to study in more 
detail by analytical means the hydrodynamic wave equa- 
tion in the axially symmetric case. We wish to find an 
explanation why at least some analytical solutions have 
been possible in this case and intend to use this insight 
to construct more solutions in a systematic way. 

In principle the collective mode problem looks very dif- 
ferent for isotropic and for axially symmetric traps: In 
the isotropic case the rotational symmetry ensures that 
angular momentum conservation gives two good quan- 



tum numbers, and therefore the wave equation is sepa- 
rable in spherical coordinates. For axial symmetry, how- 
ever, only the axial component of angular momentum re- 
mains a good quantum number, besides the energy, and 
one may expect that the system, having three degrees 
of freedom, is not integrable. In fact, this expectation is 
born out for collective excitations whose energies are nei- 
ther very large nor very small compared to the chemical 
potential (see below), which spoils all hopes to find ex- 
act analytical solutions for the modes and their spectrum 
on this energy range for axially symmetric traps. Why 
then are such solutions possible at energies in the hydro- 
dynamic regime, i.e. for energies much smaller than the 
chemical potential? 

The answer is provided by the existence, in that 
regime, of an additional conserved quantity which we ex- 
hibit explicitely below. Its existence permits the sepa- 
ration of the hydrodynamic wave equation in elliptical 
coordinates. Thereby the task of solving the wave equa- 
tion can be reduced to the purely algebraic problem of 
diagonalizing finite dimensional matrices. We use this 
method to obtain the spectrum of the low lying hydro- 
dynamic modes as a function of the ratio of the axial and 
the radial trap frequencies. 

A convenient starting point of our analysis are the lin- 
earized hydrodynamic equations as derived in Jl6[ | . They 
read 

A-kTicl , . . Ti _ . . _ 

Here tp is the phase of the macroscopic wave function, 
dp the local perturbation of the number density in the 
collective mode, po(x) is the number density of the time- 
independent condensate in the Thomas-Fermi approxi- 
mation @, p (x) = ^(/x- U(x))QQi- U(x)), where 
p is the chemical potential, M and a are the mass and 
the positive s-wave scattering length of the atoms. For 
the connection of (pi) with the Bogoliubov-equations see 
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p7[ . The trapping potential U(x) is assumed to have 
the form U{x) = (M/2)uj 2 (x 2 + y 2 ) + (M/2)w 2 z 2 . The 
surface of the condensate, in the Thomas-Fermi limit is 
defined by U = fi. Eqs. (jl|) are obtained by a gradient 
expansion and assuming Sp and to be small. Elim- 
inating from both equations and making the ansatz 
Sp(x,t) — e~ iul *iMx) one obtains the time-independent 
wave equation Muj 2 tp(x) = — V ■ (/1 — U(x))'Vip(x) 
which holds inside the condensate. We look for boundary 
conditions on the surface of the condensate which make 
the operator G — — M _1 V • (fi — J7(a:))V Hermitian, so 
that for two eigenfunctions i/j, ip with eigenvalues uj 2 , uj 2 
we have orthogonality according to 



{uj 2 - uj 1 ) / d A xip*ip = / (Fxty'Gip - i^G'f 
Jv Jv 



= — M / df{fj, - U){i)*d$/dn - ipdi>*/dn) = . 
Jav 

Here V and 91^ denote the volume of the condensate and 
its surface, respectively. Because p, — U vanishes on the 
boundary it is enough to require that ip and its normal 
derivative dip/dn remain bounded there. In the following 
it will be useful to measure lengths in units of the radial 



Thomas- Fermi radius rp = (2p/Mt 



which brings 



the wave equation into the dimensionless form uj 2 ip = Gip 
with 



G=-(uj 2 z /2)V-((l-e 2 )(l- P 2 )-z 2 )V 



(2) 



where we use cylindrical coordinates p — y 1 x 2 + y 2 , z 
and the azimuthal angle and define e 2 = 1 — uj 2 /uj 2 , 
which is positive for ujq < uj z . The operator G commutes, 
of course, with the axial angular momentum operator 
L z = —i-§^- However, there is an additional nontrivial 
operator 



B = -V 2 + {x-Vf +ix-V + e 2 d 2 /dz 2 



(3) 



which commutes with L z and with G, as one may check 
by direct calculation of the commutators. In the isotropic 
case e = 0, B may be expressed by the square of the 
angular momentum L 2 — L 2 + L 2 + L 2 and G via 

B = G/uj 2 + L 2 . Because of the existence of the three 
commuting operators G, L z , and B in the system with 
the three degrees of freedom p, z, (f) it is now manifest 
that the system is integrable, i.e. the eigenvalues of G 
can be labelled by the quantum numbers of L z and B. To 
see this explicitely we now look for variables in which the 
wave equation separates and introduce cylindrical ellipti- 
cal coordinates £, 77, which i n their oblate sph eroidal form 
JlH are defined by p = o^J (£ 2 + 1)(1 — i] 2 ) , z — cr^ry. 
These coordinates are orthogonal. Surfaces of £ = const 
are confocal ellipsoids with foci at z — 0, p — a. Sur- 
faces with constant r\ are confocal hyperboloids with the 
same foci. For uj z > ujq, i.e. < e 2 < 1, the foci at z = 
p = a are made to coincide with the foci of the ellipsoidal 



Thomas-Fermi surface if we choose a = e. The Thomas- 
Fermi surface is given by £,tf = (1/e 2 - 1) 1/2 - Then the 
interior of the condensate is described by £ in the range 
[0, £tf] an d ?7 in the range [—1,1]. 

For uj z < ujq the parameter e 2 is no longer usefull as 
it becomes negative. Instead one can define e 2 = 1 — 
uj 2 z juj 2 = —e 2 /(l — e 2 ) which lies in the range < e 2 < 
1. The foci of the T homas- Fermi ellipse now lie at z = 
±a' where a' — e'/vl — e' 2 . Therefore we now need the 
prolate spheroidal form of elliptical coordinates with 
foci at z — ±a' , p = 0. T hese coordinates are defined 
by p = & 'yj (£ 2 — 1)(1 — r/ 2 ) , z = cr'^rj where inside the 
condensate £ now has the range [1, 1/e'] while 77 has the 
same range as before. The treatments in the two cases 
are equivalent via the transformation connecting e and 
e'. In the following we shall present the equations for 
the case uj z > ujq. The final formulas for uj 2 apply for 
uj z > ujo and uj z < ujq. After the change of coordinates 
with a = e the operator G takes the form 



G = 



1 



2e 2 £ 2 + ri 2 



(4) 



(i-^(i-^))|(i-^ 2 (e 2 + i))(^ 2 + i)| 

+ (l-e 2 (e + l))l i (l-e 2 (l- V 2 ))(l- V 2 )-^ 
(l-e 2 (e + l)Kl-e 2 (l- V 2 ))(e+V 2 ), d x2 



which is now separable. The 0-dependence of its eigen- 
functions is taken care of by factors e trn ^ with the integer 
azimuthal quantum number m. Separating the operator 
in £ and r\ by making the ansatz ^>^(^)^ n (r])e' lrn ^ for its 
eigenfunctions we obtain two equations, one for 'I',, 



(i-v 2 ) T - 

dr) 



2e 2 (l-r7 2 )?7 d 



2 l-e 2 (l 

2uj 2 /uj 2 
l-e 2 (l-T] 2 ) 



di] 



-F&v (5) 



the other for ty^. It turns out that both equations are 
identical if in the equation for we substitute i£ for rj, 
i.e. ^(C) = ^j)(*0- The solution for one coordinate is 
the analytic continuation of the solution of the other from 
the real to the imaginary axis. It is easy to check that 
the separation constant (3 is just the eigenvalue of the 
operator B for the eigenfunction ^(£)$^(7/) exp(im<?!>). 
To do this one needs to express B also in the elliptic 
coordinates: 



B 



^V^{-|(i-^ + i))(^ + d| 



e + v 2 



(6) 
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The equations (|5|) contain m only quadratically. 
Therefore the energy levels are the same for ±m. Ex- 
panding tyf! in terms of associated Legendre functions 

oo 

*„ = a* p e lml (v) (7) 

t=\m\ 

we obtain from eq. (||) a second order recursion relation 
for the coefficients ag relating only even or only odd in- 
dices t. The eigenstates therefore have even and odd 
parity. The recursion relation itself is straightforward 
to obtain but lengthy and will not be written out here. 
The condition that (Q) terminates at £ max = (\ m \ + n ) 
quantizes the eigenvalue (3 

= (n + \m\)(n + \m\ +3), (8) 

which means that becomes (1 — j^ 2 )!" 1 / 2 times a poli- 
nomial of order n. 

In the isotropic case e = the operator B can be diag- 
onalizcd in spherical coordinates and its spectrum then is 
found as (3 = (2n r + £) (2n r + 1 + 3) with radial quantum 
number n r and angular quantum number I. Together 
with the connection of G and B for isotropic traps this 
gives the result of for the spectrum in the isotropic 
case. 

The eigenvalue condition for w 2 takes the form of the 
characteristic equation of a tridiagonal matrix of dimen- 
sion TV = 1 + int[n/2], which can be symmetrized by 
a suitable similarity transformation with a given diago- 
nal matrix. For fixed numbers n, \m\ we have TV dif- 
ferent solutions for w 2 , which we label by our third 
quantum number j — 0, . . . hit [n/2]. In the isotropic 
case e = the quantum number j can be expressed as 
j = mt[(£ — |m|)/2] = int[n/2] — n r as one finds by ex- 
pressing the isotropic spectrum of [FLm in terms of the 
new quantum numbers n, j < int [n/2\. 

Calculating the first levels u)(n,j,m) we get with A = 
(uj /lj z ) 2 

u 2 (0,0,to) = w 2 |m| 
w 2 (1,0,to) = oj 2 + w 2 \m\ 

lu 2 (2,j(= 0, 1), m) = Lu 2 (l + 2(|m| + 1)A 

-^^(9 - 4(|m| + 4)A + 4(|m| + 2) 2 A 2 ) 1 / 2 ) 

lu 2 (3, j(= 0, 1), m) = to 2 z ( 7 - + 2(|m| + 1)A 

-^^(25 + 4(|m| - 4)A + 4(|m| + 2) 2 A 2 ) 1/2 ) . (9) 

In the limit A -1 — > 0, which is relevant for the experi- 
ments reported in || , the mode frequencies for arbitrary 
integer n > 2j not to large can be expanded in the small 
parameter A -1 and are found as 

tu 2 (n, j, m) = io 2 (\m\ + 2j(j + \m\ + 1) + 0(A -1 )) (10) 



Remarkably, the leading order term for |m|, j not both 
vanishing is independent of n, i.e. the levels consist in 
this case of bands of closely spaced levels split only by 
small frequencies of order oj 2 /ujq. 

Eq.(||) contains, as special cases, the particular so- 
lutions previously obtained by Stringari [lq ] (the n = 
0, l;j = 0, 1 modes for all |m| and the two n = 2; j = 
0,1; m = modes). 

After our calculation had been completed and while 
this paper was being prepared a preprint became avail- 
able J2(J, in which the mode frequencies (^) where also 
found by solving the wave equation directly via a polyno- 
mial ansatz, which can be shown to be equivalent to ours. 
This method works because, as we have shown above, no 
special boundary conditions except regularity need to be 
imposed on the wave function at the surface of the con- 
densate. However, the deeper reason for the solvability 
of the equation, i.e. the additional conservation law, has 
not been identified in pofl . 

Because of the necessity to diagonalize TV-dimensional 
matrices the analytical determination of the to 2 by the 
present method for A neither very small nor very large is 
possible up to TV = 4, even though the formulas for TV = 3 
and TV = 4 are too cumbersome to be of much practical 
value. However, it is straightforward to diagonalize the 
matrices numerically and to prepare a list of the numer- 
ical values of the eigenvalues for a given ratio A, e.g. for 
the experimental value A = 1/8. In table I we give such 
a list for the 43 lowest lying eigenvalues with \m\ < 2 
in units of u>q for A = 1/8, together with their quantum 
numbers m, n, j. A numerical evaluation of the energy 
levels of the collective modes as a function of the scaling 
variable |l6| Na/do for A = 1/8 was recently reported in 
H . Here TV is the number of atoms and do = y/h/2MuiQ 
is a measure of the size of the ground state in the trap. 
The present results apply for Na/do ~ > oo. Comparing 
our results with those in [gj we find that the rate of con- 
vergence to the asymptotic case Na/do — > oo depends on 
the quantum number n and is much slower, e.g. for n = 6 
than for n = 2. For m = we find 4 additional levels 
(with quantum numbers n = 5, 7, 8, 9 and j — 0) below 
the highest m=0 level considered in ||], where the mo- 
tion of the levels was tracked from the free particle case 
to the strongly interacting case. The explanation can be 
that levels having rather high frequency in the free trap 
can move down into the considered low-energy regime as 
the parameter Na/do is increased. Following the levels 
one can expect a number of avoided level crossings due to 
the non-integrability of the system at intermediate ener- 
gies. We demonstrate the chaoticity of the corresponding 
classical system whose Hamiltonian 

H(p, x) = ^e 2 HF (p,x) - K 2 (x), 
V 2 

e„ F {p, x) = + U(x) - n + 2K{x) 

K(x) = (ji-U(x))Q(M-U(x)) (11) 
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is the Bogoliubov quasiparticle energy in local density 
approximation (Fig.l) 

In conclusion we have demonstrated that the quasipar- 
ticle dynamics in axially symmetric traps is chaotic for 
energies comparable to the chemical potential, but, for 
energies much smaller than /i we found a third conserved 
variable B commuting with the wave operator G and ax- 
ial angular momentum L z . As a consequence the wave 
equation in this regime is proven to be integrable which 
explains why some solutions have already been found in 
the literature. We have separated the wave equations in 
elliptical coordinates, reduced the eigenvalue problem to 
the diagonalization of finite dimensional matrices, and 
exhibited three integer quantum numbers m, n, j in the 
ranges — oo < m < +00, < n < 00, < j < int [n/2] 
labelling all modes and their frequencies for the experi- 
mentally realized case uj^/uiq =8. As a final remark we 
mention that we have found even the fully anisotropic 
case, where the axial angular momentum is no longer a 
good quantum number, to be integrable p2| . Two new 
conserved quantities can be identified in that case which 
reduce to L z and B in the axially symmetric limit. 
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FIG. 1. Poincare sections of the classical dynamics in cylin- 
drical coordinates p, z, (j> after the elimination of the conserved 
axial angular momentum and the azimuthal angle (f> for en- 
ergy E//J, — 1. The cut is taken at z = and diplayed in 
the scaled variables p,p p measured in units of yj2p/mu>^ and 
• s /2mp, respectively. The anisotropy is chosen as lu z /luo — 
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TABLE I. The 43 lowest levels labeled by the quantum 
numbers jmj < 2,n,j for (lu z /loo) 2 = 8 
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